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ABSTRACT 

Nonspherical collapse is modelled, under the Zeldovich approximation, by six¬ 
dimensional random walks of the initial deformation tensor field. The collapse 
boundary adopted here is a slightly-modihed version of that proposed by Chi- 
neh and Lee (2001). Not only the mass function agrees with the htting formula 
of Sheth and Tormen (1999), but the bias function and conditional mass func¬ 
tion constructed by this model are also found to agree reasonably well with the 
simulation results of Jing (1998) and Somerville et ah (2000), respectively. In 
particular, by introducing a small mass gap, we find a fitting formula for the 
conditional mass function, which works well even at small time intervals between 
parent and progenitor halos during the merging history. 
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1. Introduction 


In a hierarchical model, the gravitationally bound objects of the universe formed from 
mergers of smaller objects into larger ones. An issue of fundamental importance is to un¬ 
derstand how gravitational structures, such as galaxies and clusters, are distributed as a 
function of mass and evolving time. Press & Schechter (1974) provided a simple analytical 
model (hereafter PS model) to describe the mass distribution of dark halos (mass function). 
This theory is based on the assumption that the process of gravitational collapse can be 
approximated by spherical symmetry and it occurs when the mean density contrast of prec¬ 
ollapse regions predicted by the linear theory reaches a threshold value dc- Bond et al. (1991) 
introduced the excursion set formalism and extended the PS theory to calculations of the 
conditional mass function, which is closely related to the issue of the linear bias and halo 
merging histories. Mo & White (1996) used the extended Press-Schechter theory (hereafter 
EPS) to address the large-scale bias, which concerns how the collapsed halos trace the un¬ 
derlying dark matter. Moreover, the EPS theory is also useful for construction of the merger 
history of dark halos (Bower 1991; Lacey & Cole 1993, 1994 hereafter LC93, LC94), which 
concerns the probability that a halo of mass M 2 at time t 2 has a progenitor halo of mass Mi 
at earlier time ti. 

The PS or EPS theory has been tested against N-body simulations (Efstathiou et al. 
1988; LC94), and the comparisons suggest that the PS model overestimates the number of 
halos in low-mass regimes and underestimates it in high-mass regimes. This discrepancy 
indicates that the spherical collapse model used by the PS theory may be too simplistic. 
Many researchers have taken the nonspherical collapse process into account to correct this 
problem (e.g.. Bond & Myers 1996; Monaco 1995; Lee & Shandarin 1998; Chiueh & Lee 
2001; Sheth, Mo and Tormen 2001). Most of these works are still based on the main idea 
of the PS theory, and further incorporate either the dynamics of ellipsoids or the Zeldovich 
approximation (Zeldovich 1970). 

The main aims of this paper are to investigate the large-scale bias and conditional 
mass function by using the nonspherical collapse boundary (hereafter NCB) proposed by 
Chiueh and Lee (2001). We also test whether the results derived from nonspherical collapse 
boundary can provide better statistical predictions of dark halos than the EPS model does. 

The outline of the paper is as follows. In § 2, we summerize the PS and EPS theory, and 
briefly review the Zeldovich approximation and the nonspherical collapse boundary proposed 
by Chiueh & Lee (2001). In § 3, we show the algorithm and the results of mass function 
and large-scale bias of our model. In § 4, the conditional mass function is derived and the 
comparison with N-body simulation is presented. In § 5, the conclusion is given. 
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2. (Extended) Press-Schechter Formalism 
2.1. Mass function 


The Press-Schechter theory (1974) assumes that the initial density fluctuations are Gaus¬ 
sian random flelds. Let Sji(x) denote the overdensity held at position x in the comoving 
coordinate when the density distribution is smoothed over scale R by a window function 
Wr. The relation is 

= j 5{x')WR{x-x')(f^, ( 1 ) 

with a squared variance 

s = aHM) =< >= ^ jp(k)Wl(k)i^kHk, (2) 

where p{k) is the power spectrum of density fluctuations and WR{k) is the Fourier transform 
of the window function Wr{x — x'). The squared variance S increases as mass M decreases, 
and in this paper, we will often use S to represent the mass scale. 


Bond et ah (1991) extended the original PS method and used random walks (or called 
excursion set) to address the mass function problem. When variance S equals zero, the 
density contrast 5 is set to zero. As S increases, 5 changes randomly. When 5 first hits 
the threshold value 5c-, the random walk stops and the region of mass scale S', over which 
5 is smoothed, collapses into a bound object. This prescription is equivalent to a diffusion 
process with an absorbing boundary at <5 = 5c-, which has been solved by Chandrasekhar 
(1943) when the window function is chosen to be a sharp k-space Alter. Thus, one can 
calculate the probability of a trajectory having its first hit at the threshold 5c in the mass 
range between S and S' -|- dS* to represent the mass function, which is described by 

where the functional form f{x) represents the differential probability in the range between 
X and X -|- dx- Moreover, the (comoving) number density of halos with mass M is 


dn 

dM 


dM 


^f{S,5c)\^\dM, 


( 4 ) 


where pf, is the mean mass density of the universe. Alternative choices of window functions, 
such as real-space top-hat Alter and Gaussian Alter, are also often used to determine the 
relation of the squared variance S and mass M. The localized Alters yield correlated steps in 
the above random-walk framework. Comparisons among different types of window function 
are well discussed by Bond et al. (1991), LC93 and LC94. 
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For simplicity, we define a variable v = 5c/8^^“^ and one can easily show that the mass 
function can be represented as 


ijf{ij)du = 




(5) 


Here 5c is often chosen as 1.69, which is predicted by the spherical collapse model. The 
main advantage of the excusion set approach is that while the mass function depends on the 
variance S and the value of collapse threshold 5c-, the cosmological models enter this problem 
only through the relation between S and the mass M. 


Upon comparing with N-body simulations, one hnds that the EPS formalism overesti¬ 
mates low-mass and underestimates high-mass halos. Of course one may vary the value of 
the collapse threshold 5c-, but it is not possible to match simulation results well for all mass 
ranges. 


2.2. Conditional mass function 


In the picture of the excursion set theory, the conditional probability represents the 
fraction of trajectories hrst hitting a higher threshold 5ci among those trajectories that start 
from a lower threshold 5c2- In the framework of EPS theory, Lacey & Cole (1993) proposed 
that the conditional probability can assume a modihed form where the origin of the mass 
function (0,0) shifts to (5*2, 5c2)- Therefore, the conditional probability is simply 


f{S,,5ci\S2,5c2)dS^ = 


((5cl — 5 c 2 ) 


V^(^1 - 

Here, it demands > S 2 and 5ci > 5c2- 


exp[ 


(5cl - 5 c2)' 

2(^1 - S2) 


■]dSi- 


( 6 ) 


Conditional probability is useful for investigation of large-scale bias as well as merger 
events of halos. In the hrst case, the physical meaning of 5c2 is not the collapse threshold 
at some specihc time, but denotes the mean density contrast of a large uncollapsed region 
of mass scale S 2 - Thus we shall change 5c2 into ^2 to avoid confusion. If we dehne the bias 
parameter b to be the fraction of excess dark halos with mass scale Si that can be found 
in a larger mass scale S 2 , then under the limit S 2 ^ Si and 52 -C 5ci, Mo & White (1996, 
hereafter MW96) use the EPS formalism to derive that 

f - 1 vl-1 

5ci 5ci 



where vi = 5ci/S\^‘^- For convenience, we dehne a characteristic mass M*, which satishes 
the requirement S{M*) = 5/^- When Si < 5/^, i.e.. Mi > M*, the bias factor b is positive. 
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which means the density of halos with masses greater than M* is larger in the denser region 
of dark matter. On the other hand, when Si > or Mi < M*, the opposite holds. 

In the case of merger events, 6 c 2 denotes the collapse threshold at a later time t 2 (as 
compared with ti). Again, the meanings of notations change: Sd —^ Sc{zi) and Sc 2 —^ Sc{z 2 ), 
the collapse thresholds at redshifts zi and Z 2 , respectively. Thus, f{Mi,Sci\M 2 ,Sc 2 )dMi 
represents the mass fraction of M 2 at redshift Z 2 that is in the form of collapsed halos of 
mass Ml in an earlier epoch Zi. 


2.3. Zeldovich approximation 

Zeldovich (1970) provided a good way to describe the development of density perturba¬ 
tions into the nonlinear regime. Let a{t) denote the uniform expansion of the background. 
The physical coordinate r and the comoving coordinate x are simply connected by 


r = a(t)x. 


( 8 ) 


In the Zeldovich aprroximation, the motion of a particle is described by 


x = g + b(t)p{^, 


(9) 


where q is the Lagrangian coordinate (initial location) of a particle. The second term is 
the perturbation of a particle in the Lagrangian coordinate. From the conservation of mass 
condition, we have 

p(T, t)(fx = p(g, t)(fq = p(g, Q)Sq. ( 10 ) 


It implies 


p(f, t) 


_dSq _ 1 _ _ 1 

^ d^x ^ det{^) ^ detj' 


( 11 ) 


The tensor J in the denominator is the sum of an identity tensor and a deformation tensor 
b. Let p(^ = —V^(^,then 


Dij = hit) 


dpi 

dqi 


-hit) 


dqidqj 




( 12 ) 


Because D is real and symmetric, it has three real eigenvalues, A = (—Ai, —A 2 , —A 3 ), at 
every spatial location. Thus, Eq. (11) can be rewritten as 

p{x, t) = P' 


(1 — Ai)(l — A 2)(1 — A 3 ) 


(13) 



- 6 - 


When A's are small, the density contrast 6 is 

<5 = ^ (Ai + A2 + Ag). (14) 

P 

The sign of the value of A determines the divergence (negative A) or convergence (positive 
A) of the local motion. Some examples illustrate the physical meaning of Eq. (13). If 
Ai = A 2 = Ag, when all A’s approach unity the density becomes inhnity, and it implies the 
spherical collapse case. If a certain A* grows to unity before the other two A’s do, then the 
local ellipsoid will collapse into a pancake. It is easy to realize that the latter be much more 
probable than the former, and therefore the gravitational collapse should be nonspherical. 

Since the deformation tensor D is symmetric, it has 6 degrees of freedom: du, ^ 22 , dgg, 
di 2 , d 23 and dig. The first three represent the diagonal terms, and the rest represent the 
off-diagonal ones. In addition, one can construct 6 independent Gaussian variables, namely, 
^ = (di) 1/2, Us, Vii 1/5, y%)- The 6 components of the deformation tensor can be interpreted as 
linear combinations of yi ^ ye'. 


3 1 

dll = (yi + ^y 2 + yg^3)/3. 

(15a) 

2 

d 22 = -{yi - ^|/g)/3. 

(15b) 

3 1 

d33= (yi ^^2+^d3)/3, 

(15c) 

, 1 
di2 — -^yi) 

15 

(15d) 

d23 = Tyds, 

15 

(15e) 

di3 = Tyd6- 
15 

(15f) 


2.4. Nonspherical collapse boundary (NCB) 

In order to resolve the discrepancy of the mass functions between the PS prediction and 
the N-body simulation, Chiueh & Lee (2001) investigate how the nonspherical collapse pro¬ 
cess may be relevant. In addition to adopting the fundamental concepts of the PS hypothesis 
and the excursion set approach, they have noticed that the collapse condition should contain 
full information of the three eigenvalues (Ai, A2, Ag), and that the nonspherical effect is more 
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dominant in the low-mass regime than in the high-mass regime. Taking these properties 
into account, they propose a nonspherical collapse boundary (NCB hereafter) having the 
following form: 

;f = (l + ^)'’. ( 16 ) 

Uc P 

where /3 is a htting parameter taken to be 0.15, 5c remains the collapse threshold predicted 
by the spherical model, which is taken to be 1.5 by Chiueh & Lee (2001) rather than 1.69 
due to the correction given by Shapiro et ah (1999), and the dehnition of is 

r‘ = l[(Ai - \i? + (Aj - As)" + (As - A,)"]. ( 17 ) 

In the case of spherical collapse, the three eigenvalues equal to each other. Therefore, the 
quantity r can be regarded as an indicator for the degree of nonsphericality. 

The form of NCB has several advantages. First, this boundary is smooth and rotational 
invariant in the eigenvalue (Ai,A 2 ,A 3 ) space. Second, when r is close to zero, the NCB is 
reduced to the spherical collapse condition as expected. Third, the variable is related 
to the angular momentum (Catelan & Theuns 1996), so that the NCB provides additional 
information about the nonspherical collapse halos. 


3. Mass Function and Bias for NCB 
3.1. Mass function 

Here we slightly improve the nonspherical collapse boundary of Chiueh and Lee (2001) 
into this form: 

where 5c = 1.5654, Pi = 0.26 and P 2 = 0.159. The three symbols 5c, Pi and P 2 are all htting 
parameters to obtain a mass function closest to the result of N-body simulation (Sheth & 
Tormen 1999, hereafter ST99). Fine tuning of these parameters would not change the entire 
shape of mass function too much but it is necessary for us to perform the bias problem. The 
eigenvalues A satisfy the eigenvalue equation 

A^ -h (dll + *^22 + -\- {diid22 + <^22^33 -|- d33dii — d^2 ~ *^23 ~ *^ 13 )-^ ~ 

(dlld23 -|- (^22^13 -|- d33di2 — ^11(722^33 — 2(712^23^13) = 0. (19) 


Therefore, 


Ai + A2 -l- A3 — —{dll + d 22 + *733), 


(20a) 



(20b) 

(20c) 
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A1A2 + A2A3 + A3A1 = {diid ‘22 + <^ 22 C ?33 + <^ 33<^11 ~ <^12 ~ ~ *^13)) 

A 1 A 2 A 3 = ((^11(^23 + '^22'^13 + <^33'^12 ~ '^ll'^22<^33 ~ 2(ii2(i23<^13)• 

Using Eqs. (15), (20a) and (20b), it is easy to convert into a function of |/i, I/ 2 , 1 / 3 , 1 / 4 , 1 / 5 , 
and y^. The result is 

= ^{yl + yl + yl + yl + vD- ( 21 ) 

This neat equation has been pointed out by Sheth & Tormen (2002). Combined with 5 = yi 
(from Eqs. (14), (15) and (20a)), we can simulate the random process in the 6 -dimensional 
U-space with an absorbing boundary, instead of conducting the equivalent, but more com¬ 
plicated, random walks in (Ai, A 2 , A 3 )-space (Chiueh & Lee, 2001). 

In Fig. 1 we compare the mass function derived from the U-space diffusion with two 
other types of mass functions. The hlled circles are the data of our six-dimensional random 
walks, the dashed line shows the EPS’s analytical prediction, and the solid line shows the 
Sheth & Tormen’s mass function (ST99), which has the form 

^/(O = 2A11 + (22) 

where a = 0.707, q = 0.3 and A = 0.322. We conclude that our results agree quite well with 
that of N-body simulation given by ST99 (c.f., Eq. (22)), except at very high mass regime. 
There are two reasons for us to retain this discrepancy at high mass. First, the statistics 
of very high mass halos is not well-controlled in N-body simulation. Therefore one can not 
reliably specify the number of such halos. Second, it is believed that the very high mass 
halos should follow spherical collapse process which is closer to the EPS method, and our 
results show this tendency. As a further comparison, we also demonstrate in Fig. 1 the mass 
function given by the original nonspherical collapse boundary proposed by Chiueh & Lee 
(2001). It is shown that the two NCB mass functions are roughly the same except that our 
mass function is more accurate in the intermediate mass range. Just as we have mentioned 
earlier, the high accuracy is necessary when one wants to derive the bias factor precisely 
since the bias factor is related to the small difference between the biased mass function and 
the unbiased mass function, an issue to be addressed in the next section. 


3.2. Bias 

The peak biasing was originally introduced by Kaiser (1984) to explain the clustering 
of Abell clusters, and the same idea has also been used to discuss whether the galaxy dis¬ 
tribution traces the mass distribution (Bardeen et ah 1986). Despite focusing on specihc 
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bound objects, the bias parameter also reveals the spatial clustering of dark matter halos. 
One common dehnition of the bias in the Lagrangian coordinate is 


6n ^6p 
n p ’ 


(23) 


or 


^ n(M, = (50 - n(M, = 0) 

^ ^ n{M,6,\6 = 0) X S' 


(24) 


where n{M, (5c 1(5) represents the number of bound objects with mass M which can be found 
in the large-scale uncollapsed region of density contrast S. The second term of the numerator 
on the right hand side n{M,Sc\S = 0), or named ’’the unbiased mass function”, is actually 
the conventional mass function we have been discussing so far, and n{M,Sc\S = S') is ’’the 
biased mass function” in the sense that it is embedded in the uncollapsed region with hnite 
density contrast S'. Mo & White (1996) have derived the bias factor under the limit S « Sc, 
as shown in Eq. (7), and they also give another equivalent expression for the bias evaluated 
at a sufficiently large separation r: 




(25) 


where ^(r) is the two-point correlation function and the subscripts hh, mm denote the halo- 
halo and matter-matter components respectively. This dehnition is often used to hnd the 
bias factor in N-body simulation. 

In the EPS picture, one can calculate b as: 

_ n(M, Sc - (5'|(5 = 0) - n{M, (5c|(5 = 0) _ 1 dn{v) _ 1 d(lnn(z/)) 

n{M,Sc\S = ^) X S' a{M) X n{v) dv Sc dilnv) ’ 

indicating that the bias function is directly related to the derivative of mass function in 
the log-space. However, this simplicity is no longer valid in the NCB approach due to the 
additional r-distribution, and therefore the bias function can be another test for the goodness 
of the NCB. The biased mass function n{M, (5c|(5 = S') can still be obtained from the excursion 
set approach. The natural choice for the initial condition of six-dimensional random walks 
for the biased collapse is that yi = S', whereas 1/2 to i/q are Gaussian distributed random 
numbers with the same variance equal to S'. Comparing the biased mass function with the 
unbiased mass function already obtained in § 3.1, we are able to compute the bias function 
dehned in Eq. (24). 

In Fig. 2, we show the bias functions derived from our approach for three different 
initial power spectra of indices, n = 0, —1, and —2. The hlled circles show the results of our 
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algorithm of Y -space random walks. The solid curves denote the fitting formula of N-body 
simulation results found by Jing (1998): 

h{M) = (^ + l)( 0 - 06 - 0 . 02 n)(^ ^ (27) 

y Oq 

As a comparison, Fig. 2 also shows the predictions of the EPS model (MW96) and the 
moving-barrier model (Sheth, Mo & Tormen 2001). Although there exit statistical fluctua¬ 
tions in our data, the filled circles show that the bias function generated from the nonspherical 
collapse boundary is quite reasonable and has a similar tendency as the results derived from 
the moving barrier model in the low-mass scale. But the curves of our results turn to coincide 
with the EPS prediction at high masses. This tendency is also seen in Jing’s simulations. 


4. Conditional Mass Function 

According to the hierarchical scenario, dark matter halos are formed through merging 
with comparable-size halos and accretion of small halos. All kinds of interactions between 
halos result in changes of gravitational potential and consequently the behavior of gas com¬ 
ponents. Therefore, merger/accretion rates are believed to play a crucial role in the galaxy 
formation. Observations show that merger rates increase with redshift (Carlberg et al. 1994; 
Yee & Ellingson 1995), despite the different definitions of the merger rate. On the other 
hand, the dependence of merger rates on redshifts and on environments has also been stud¬ 
ied with N-body simulations (e.g., Gottlober et al. 2001), and the results roughly agree 
with the observations. It is expected that the semi-analytical approach (e.g., EPS model 
and NOB model) can offer succint description for merger histories. However, we must keep 
in mind that the semi-analytical methods deal with only the isolated halos, discarding the 
information of substructures within each halo. The attempt of explaining galaxy formation 
inside galaxy groups or clusters with semi-analytical approaches must be modified. Despite 
such a limitation for the cloud-in-cloud problem, the merger-tree technique is often needed to 
establish merger histories in the field ( Kauffmann & White 1993; Somerville & Kolatt 1999). 
The merger-tree method is a Monte-Carlo method based on the conditional mass function. 
In the literature, EPS’s conditional mass function is often used to construct merger trees. 
However, the results indicate that there still exhibit discrepancies with the merger history 
extracted from N-body simulations, which is thought to be caused mainly by the EPS model 
rather than the specific features of merger trees (Somerville et al. 2000). Therefore, in 
this section, we use the NCB model to generate a new conditional mass function with the 
six-dimensional random walk procedure and derive a fitting formula for it. 
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4.1. Construction of conditional mass function 


There are two scenarios, the active picture and the passive picture, to represent the 
diffusion process for epochs earlier than the present. The former is with a fixed collapse 
boundary and let the variance of some particular mass grow with time. The latter is just 
the opposite which fixes the variance of mass M and lowers the collapse boundary with 
increasing time. In order to simulate the merging problem, we adopt the second picture for 
simplicity. For a given redshift 2 :, the collapse boundary condition can be rewritten as 


5c{z) 


(1 + 


1 r 


= 0 . 26,/32 


0.159) 


(28) 


where Sc{z) denotes the vertex of the boundary at redshift 2 : and Sc{z = 0) is just Sc given 
in Eq. (18). The evolution of Sc{z) is inversely proportional to the growth factor of the 
variance a. For example, Sc{z) oc (1 + ;^) in S'CDM model. As shown in the Fig. 3, boundary 
1 corresponds to the collapse condition at a time earlier than boundary 2. That is, Eq. (28) 
is nothing but to say that boundary 1 should be stretched both in r and S axes relative 
to boundary 2 to preserve the self-similar evolution of fiuctuations. After collecting the 
trajectories that start from the origin and collapse at boundary 2 , we continue the random 
walks until they hit boundaryl and are absorbed. 


Let f{Si, (5 c(2:i)|>S'2 , Sc{z 2 )) denote the fraction of trajectories hitting boundary 1 at mass 
scale Si to the total trajectories hitting boundary 2 at mass scale 82 - The physical meaning 
of /(S*!, (5 c(2:i)|*S'2 , (5 c(2 :2)) is the probability that a point in a halo with mass scale S 2 at 
redshift Z 2 was in a halo with mass scale Si at redshift zi. As mentioned earlier, although 
S is the variance of the density fluctuation corresponding to a specific mass range, we use 
it to interpret the mass scale over which the density contrast S is smoothed. Apparently , 
Si is always larger than S 2 to ensure that the mass of halo increase with time. Moreover, 
due to the self-similar scaling property of the nonspherical collapse boundary, / does not 
depend on time explicitly, but only on the ratio of the separation, namely Sc{zi)/Sc{z 2 ) = n. 
The actual time evolution of the conditional mass function is hidden in the variable fi or fj,' 
defined below, which is a combination of mass and time. Therefore, by simulating random 
walks with a particular chosen Z 2 = 0, we are able to construct the conditional mass function 
which can also be applied to any other Z 2 - 

Recalling the Pess-Schechter’s boundary, its conditional probability is given in Eq. (6), 
or in a concise form, 

= 2 X (^)^/2exp(-^)d/i, (29) 

27r 2 

where /i = { 6 c{zi) — Sc{z 2 ))/{Si — S' 2 )^/^ = ASc/ (AS*)^/^. Apparently, for given zi and Z 2 , the 
conditional probability is a function of only the difference of the variance S in EPS picture. 
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However, the conditional probability actually depends also on the initial mass scale S 2 in 
our scenario. This is due to that the separation of two boundaries is S' 2 -dependent, as is 
obviously indicated in Fig. 3. 

Here we simulate several conditions for different boundary separations: k = 1.1, 1.3, 
2.0, 3.0 and 4.0 with Z 2 = 0. It helps the construction of a new conditional mass function by 
considering the following two extreme cases: when k approaches unit (the two boundaries are 
close), the conditional mass function is expected to be close to the extend Press-Schechter’s 
case because the two boundaries are almost parallel and most points collapsing at the first 
boundary will soon collapse in few steps at the second boundary. On the other hand, when k 
is very large, points in the second boundary relative to the first boundary are like that they 
are located at the origin, so that the conditional mass function should behave like the Sheth 
and Tormen’s mass function (hereafter ST mass function). We hnd that our simulation data 
can still be fitted by the following function: 

.. /2 /2 

= 2A{k){1 + —)(^)^'^^exp(-^)dp', (30) 

jj, 2% 2 


where 


' _ [Sc{zi) - Sc{z2)]e{S2, k) 

(^1-^2)1/2 

0.178 


A = 0.322 + 


(31) 

(32) 

e(+ 2 , k) = e(x) = 0.036a;^ - 0.309a;3 + - 1.060a; + 1 (33) 

with X = {\fS 2 / 5c{z2) — 0.25)//s;, and qi^n) is required to satisfy the normalization condition 


K 


A = 



(34) 


Further explanation of the htting procedure for these htting parameters is described in 
Appendix A. 


The mathematical form of Eq. (30) is similar to the ST mass function, which demands 
A, q to be constants. Here we extend the ST mass function to the conditional mass function 
and introduce A{k), q{K) and e(S' 2 , k) as fitting functions. The dehnition of fj! here has been 
so motivated to £t the probability distribution derived from random walks. 


However, it has been noted that when the time interval is small, i.e., Z 2 close to Zi, 
the concept of random walks fails to describe the merger events because the excursion set 
neglects the correlations between scales (Peacock & Heavens 1990, Sheth & Tormen 2002), 
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and this explains why the conditional mass fnnction derived from the moving-barrier in a 
recent report (Sheth & Tormen 2002) may not be applied to constrnct the merger tree. 
Althongh Eqs. (30) and (31) faithfnlly describe the resnlt of onr random-walk processes, in 
order to make the conditional mass fnnction usefnl even at small time interval, we hnd that 
fj,' shonld be modihed into 

/ ^ [5c{zi)-5c{z2)\e{S2,K) 

^ ~ {Sr-S2-v5l{z2)Y/^ ■ ^ ^ 

The replacement of S'i(Mi) — S' 2 (M 2 ) by 5'i(Mi) — S' 2 (M 2 ) —r]5l{z2) means that there exists 
a small mass gap between parent halos and progenitor halos. This small mass gap only 
becomes important when the look-back time is small and becomes negligible for a large 
look-back time. A small positive rj serves to correlate the probability of progenitor halos 
to that of parent halos. We hnd that rj = 0.06 gives the best result when compared with 
simulations, as discussed below. Furthermore, the variance and threshold value appear as 
'/S/5c{z2) in fi', such that the active picture and the passive picture are consistent with each 
other. 


4.2. Comparison with N-body simulation 


In order to test whether the conditional mass function derived in the previous sec¬ 
tion is valid for describing the formation history of dark halos, we compare our conditional 
mass function with the results of the particle-particle/ particle-mesh simulation given by 
Somerville et ah (2000). This simulation is a run for the rCDM model (hi = 1, T = 0.21, h = 
0.5) with box size 170 Mpc, particle number 256^. The halos are identihed by using the 
’friends-of-friends’ (FOF) algorithm. 

To convert the conditional mass function into the probability that a halo with mass M 2 
at redshift Z 2 had a progenitor in the mass range Mi to Mi + dMi at redshift Zi, we have 
used (LC93) 


dn{Mi zi\M 2 ,Z 2 ) ^ Mlf^Si,5Yzi)\S2,SYz2))\^\dMi, 

dMi m / ^ ^ ^ ’’^dMY 


(36) 


where 

f(Si,S,(zi)\S,,h(z2)) = (37) 

The mass range of the parent halo chosen by Somerville et ah (2000) is 7.9 ~ 12.OM^ at the 
present time {z 2 = 0), where Ml ~ 5.0 x IO^^Mq. Because the parent mass range is wide, we 
adopt a fair assumption that the parents of different mass obey the ST mass function. We 
use this distribution of M 2 to weigh the conditional mass function for the parent population. 
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In Fig. 4, we compare the conditional mass functions predicted by our nonspherical collapse 
procedure and by EPS with the results of the rCDM simulation (Somerville et ah, 2000). 
The four panels correspond to 4 different redshifts: Zi = 0.2, 0.5,1.1 and 2.1. The results of 
rCDM simulation are represented by hlled triangles; the dotted curves denote predictions of 
the EPS theory; the dashed curves show predictions by using Eq. (30) with the uncorrected 
dehnition of p' (c.f., Eq. (31)). The solid curves, which are closer to the simulation results 
are our predictions with the corrected dehnition of fi' (c.f., Eq. (35)). It can be seen that 
introducing the small correction r] in Eq. (35) helps reproduce the conditional mass function 
more accurately for massive progenitor halos at small zi. The solid curves and dashed 
curves can almost coincide when zi gets larger, as the effect of r] becomes negligible. The 
comparison shows that our conditional mass function (c.f., Eq. (30) & Eq. (35)) provides a 
better description than EPS does for a wide range of redshift. The successful description of 
our conditional mass function suggests that the nonspherical collapse boundary can be more 
’’realistic” and may be applied to construction of merger trees. 


5. Conclusion 

In this paper, we have attempted to show that the nonspherical collapse boundary 
(NCB) reproduces the mass function, bias factor and conditional mass function better than 
the Press-Schechter theory. This model successfully addresses these subjects over a wide 
range of mass scales, while the (extend) Press-Schechter theory always fails in some mass 
range. 

Through the simulations of 6-dimensional random walks, we have shown that the bias 
function generated by NCB agrees with the N-body simulations reasonably well. Because 
NCB is more complicated than the PS’s collapse boundary, the relation between the bias 
function and the conditional mass function is no longer as simple as the EPS model. Espe¬ 
cially, the two issues are not of identical representation in the framework of NCB model. We 
provide a formula for the conditional mass function of the nonspherical collapse boundary, 
which can be applied to any two redshifts. In order to overcome the problem that random- 
walks fail to describe the merger history at small look-back time, we add a minor correction 
in the definition of variable p' (c.f. Eq. (35)) while maintaining the overall mathematical 
form intact. This formula can describe the conditional mass probability more accurately 
than the EPS formalism over a wide mass range, even when the two epochs are close. All 
the results we obtained from NCB suggest that NCB is sufficiently reasonable to reproduce 
the statistical properties of dark halos. 

There are two interesting issues pertinent to NCB; one is its natural association with 
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the halo angular momentum through its dependence on the quantity r. The halo angular 
momentum distribution predicted by the NCB model is reported elsewhere (Chiueh, Lee 
and Lin, in preparation). Second, as mentioned in § 4, the conditional mass function offers 
an useful tool to construct the merger history through merger trees. Recent works have 
indicated that the merger history constructed by using EPS’s conditional mass function 
deviates from that extracted from N-body simulations, and it is thought to be caused mainly 
by the spherical-symmetry assumption of the EPS model rather than the specific features 
of the merger tree scheme. Thus, we expect the agreement can be greatly improved by the 
NCB conditional mass function given in this work. 


We thank Dr. Y. Jing for useful comments on this paper. We also thank Dr. H. Mo 
and Dr. R. Wechsler for pointing out typos. This work is supported in part by the National 
Science Council of Taiwan under the grant, NSC-90-2112-M-002-026. 


A. The fitting procedure of conditional mass function 


After performing the random walks between two collapse boundaries, which are men¬ 
tioned in § 4.1, we use the mathematical form of ST mass function and tune the values of 
A, q and £ to adjust our results of random walks. It should be emphasized that A and q are 
not tuned independently because the integral of the function should be normalized to 
unity, due to the assumption that all the trajectories starting from the first boundary will 
eventually hit the second boundary. Thus, we have 

/2 

J = j 2 A(k)(1 + = 1, (Al) 

then we get 



In addition, it is expected that all three parameters A, q and e should yield a fitting formula 
of conditional probability to match the two extreme cases, as mentioned in § 4.1. For k 
approaching unit, i.e, the two boundaries are very close, A{k) gives the value and q is 
thus zero, which agrees with the Press-Schechter’s conditional probability. For k to be very 
large, A and q approach 0.322 and 0.3 respectively, the values given by Sheth & Tormen 
(1999). The results of our 6-dimensional simulation and our fitting formula Eq. (30) and 
(31) are shown in Fig. 5 ~ Fig. 9. The five figures correspond to different separations of 
boundaries, i.e., k = 1.1,1.3, 2.0, 3.0, and 4.0. The top panel of each figure is the result for 
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initial mass scale = 0.7, the middle panel is for ^/S 2 = 1.5 (around M*) and the bottom 
is for = 2.9. The opaque circles represent the data from 6-dimensional simulation and 
the solid lines are the best fit using Eq. (30) and (31). 

We hnd that for the case k is hxed, A does not vary with the initial mass scale S 2 , 
i.e., A is a function of k only. In Fig. 10, we show the relation between A and 1/k. The 
square symbols show the htting values of A for five different k plus two boundary values: 
{1/k, A) = (0,0.322) and (1,0.5), which recover the limit conditions associated with EPS 
and ST cases. The solid line represents here is a linear relation between A and 1/k\ 

A = 0.322+^^^, (A3) 

K 

which is determined by the above two boundary points. This line almost passes through 
the other five htting points and thus is a good representation for A{k). Combined with Eq. 
(A2), q can be determined immediately for each k. 

Our tests also show that e is a function of a combined variable x = {^/^/5c{z2)—13■23)/K. 
The best ht of the function is (see Fig. 11) 

e{x) = 0.036a;^ - 0.309a;^ + OMAx^ - 1.060a; + 1. (A4) 

This result is consistent with the geometry of two boundaries. Since we have the r-distribution 
of various mass scales S 2 A that hit the hrst boundary, one can estimate roughly the distance 
between two boundaries for each mass scale. It is easily seen that for a hxed k, the gap d 
initially decreases to reach a minimum and then rises when S 2 continues to increase. The 
behavior of d/{Sc{zi) — Sc{z 2 )) agrees with that of e dehned above. 
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Fig. 1.— Comparison of the mass function derived from EPS model (dashed line) and the 
nonspherical collapse boundary (hlled circles) with Sheth & Tormen’s htting formula of N- 
body simulations (solid line). The opaque triangles show the mass function generated from 
the original nonspherical collapse boundary proposed by Chiueh & Lee (2001). 
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M/M* 



M/M* 

Fig. 2.— The bias factor as a function of M/M* for three different initial power indices:(a) 
n=0, (b) n=-l and (c) n=-2. The circle symbols show the results of our algorithm of F-space 
random walks. The solid curves are Jing’s fitting formula derived from N-body simulations. 
The long dashed lines denote the analytical bias prediction based on EPS (MW96), and the 
dot-dashed lines show the results predicted by moving-barrier model (Sheth, Mo & Tormen 
2001 ). 
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Fig. 3.— The nonspherical collapse boundaries of two different epochs. The higher the 
boundary is, the earlier the collapse time. 
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Fig. 4.— Comparison of the conditional mass function predicted by EPS (dotted curves), 
uncorrected NCB (dashed curves) and corrected NCB (solid curves) models with the results 
of the rCDM simulations (triangle symbols). The parent halos are chosen at z = 0, and the 
progenitors are at (a) 2 : = 0.2, (b) 2 : = 0.5, (c) 2 : = 1.1 and (d) 2 : = 2.1. 
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Fig. 5.— The conditionalmass function for k, = 1.1. The opaque circles represent the results 
of the six-dimensional random walks and the solid lines are the best £t using Eqs. (30) and 
(31). The mass scale at time Z 2 is chosen as (a) = 0.7, (b) = 1-5 and (c) = 2.9. 
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Same as the previous figure, but for k = 1.3. 
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Same as the previous figure, but for k = 2.0. 




(S,-S,)xf(S,-S,) 





(S,-S,)xf(S,-S,) 



Same as the previous figure, but for k = 4.0. 
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A 

0.5 



Fig. 10.— A — 1/k relation. The square symbols show the hve htting values of A for our 
data, the hlled triangle represents the EPS case, the hlled circle represents the case of ST 
mass function, and the solid line denotes the linear relation between A and 1/k (c.f., Eq. 
(32)), 
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Fig. 11.— e factor as a function of x{= 5c{z2) — 0.25)/fi;). Similar to Fig. 10, the 

square symbols show the htting values of e for our data and the solid line denotes the htting 
relation of the square symbols (c.f., Eq. (33)). 




